// main_fy3dac_mosiac.cpp
// 输入不同坐标的输入数据，拼接输出等经纬数据。
// 2024/4/17
// 
// 这里假设数据经过了云和阴影掩膜过（mask_out_cloudNshadows），留下的都是好数据，
// 只做最小蓝光合成，不考虑百分比，最小第几个什么的。
// 2024/6/24
// ------------

#include <vector>
#include <string>
#include "wdirtools.h" //../../sharedcodes/wftools.h"
#include "wstringutils.h"
#define WGDALRASTER_H_IMPLEMENTATION
#include "wGdalRaster.h"
#include <iostream>
#include "gdal_priv.h"
#include <cmath>
#include <cassert>
#include "filenameutil.h"
#include "getopt_pp.h"
#include <unordered_map>
using namespace std;
using namespace GetOpt;



int main(int argc , char** argv )
{
	cout<<"A program to mosiac multi-bands input data into one WGS84 output data. 2024/4/17"<<endl;
	cout<<  "version v1.0.0.0"<<endl ;
	cout<<"version v2.0.1.1 rename to minblue_mosiac 2024-4-24"<<endl;
	cout<<"version v3.0.0.3d add extra band for int tag; add cloud shadow remove."<<endl;
	cout<<"version v3.0.1.0d use neighbor pixels for shadow testing."<<endl;
	{
		// 我判断云阴影的原理：假设像素A在8天过程中晴天和有云的概率各是50%，那么云阴影的概率肯定极小，
		// 云阴影的成立条件是边上像素是云（0.5概率），另一侧是晴空（0.5概率），那么A的云阴影的概率应该是 0.5*0.5=0.25的概率。
		// 因此我这里 v3.0.0.x 这个版本采用云阴影概率最大不超过 25%计算。比如8天的话，最多2天是阴影，那么我就可以取最小的三个值，从中选择一个不是阴影的最小值。
		// 这里我的算法记录如下：设总天数N，阴影最大天数S，最小值采样个数MI。
		// 如果N<=5,那么S=1，MI=2；
		// 5<N<=10,那么S=2，MI=3；
		// N>10 那么S=N*0.2， MI=S+1；
		// 关于最后结果的选择，现在只考虑MI=3，如果MI-1，MI-2，MI-3相近，那么取值MI-1
		// 如果MI-2和MI-3相近，那么取值MI-2
		// 如果MI-1和MI-2相近，那么取值MI-3
		// 上面这个套路太复杂了，简化如下，保留最小blue的两天数据MI-A和MI-B，如果A和B相差在+/- 15%以内，那么就取偏小那个值，反之取偏大那个值。
		// 记录一下v3.0.0.3d 问题如下，两个最小值不能排除阴影，基本把两个最小值的差异缩小到20以内也不行。
		// 所以为了改进我打算在判断云上下功夫，比如确定某个像素是云，那么我干脆排除掉云周边 5 个像素缓存区的像素。 (这个不行2024/4/29 13:56)
		// v3.0.1.x 这里计划这么判断阴影，逐个像素判断，到了像素A，判断他太阳方向上是否有云，如果有云那么直接放弃这个像素。
	}
	cout<<"version v4.0d use cloudshadows mask file to remove clouds and shadows pixels. 2024/5/31"<<endl;
	// v4.x give up v3.x algorithm. Just use clouds(1) and shadows(2) mask files.
	// a input file must has a cs file, otherwise the input file will be get rid of.
	cout<<"version v5.0.0.1 we believe the input data are already mask out clouds and shadows,"<<endl;
	cout<<" this program only do minblue compose. 2024/6/24"<<endl;

	cout<<"usage: minblue_mosiac "<<endl;
	cout<<"    -i input dir "<<endl;
	cout<<"    --prefix filename prefix"<<endl;
	cout<<"    --tail filename tail"<<endl;
	cout<<"    --datepos date str pos in filename"<<endl;
	cout<<"    --datelen date str len in filename"<<endl;
	cout<<"    --datemin min date included"<<endl;
	cout<<"    --datemax max date included"<<endl;
	cout<<"    -b bluebandindex0based "<<endl;
	cout<<"    --vmin validbluemin --vmax validbluemax "<<endl;
	cout<<"    -f fillvalue "<<endl;
	cout<<"    --oullr out extent in ulx uly lrx lry"<<endl;
	cout<<"    -o output filename. output have same bands, srs and resolution with first input."<<endl;

	GDALAllRegister() ;

	GetOpt_pp ops(argc,argv) ;
	string indir;
	string fnprefix, fntail ;
	int dpos(-1) , dlen(0) , dateMin(0) , dateMax(0) ;
	int blueIndex0B = 0 ;
	double blueMin = 0 ;
	double blueMax = 0 ;
	double fillvalue = 0 ;
	string outfilename ;
	vector<double> outExtent_ullr ;
	ops>>Option('i',indir)>>
	Option("prefix",fnprefix)>>
	Option("tail" , fntail )>>
	Option("datepos",dpos)>>
	Option("datelen",dlen)>>
	Option("datemin", dateMin)>>
	Option("datemax", dateMax)>>
	Option('b', blueIndex0B)>>
	Option("vmin", blueMin)>>
	Option("vmax", blueMax)>>
	Option('f',fillvalue)>>
	Option("oullr", outExtent_ullr)>>
	Option('o',outfilename);

	if( indir=="" ) {
		cout<<"Error. empty indir"<<endl;
		return 11 ;
	}
	if( fnprefix=="" ) {
		cout<<"Error. empty fnprefix"<<endl;
		return 11 ;
	}
	if( fntail=="" ) {
		cout<<"Error. empty fntail"<<endl;
		return 11 ;
	}
	if( dpos<0 ) {
		cout<<"Error. empty datepos"<<endl;
		return 11 ;
	}
	if( dlen<=0 ) {
		cout<<"Error. empty datelen"<<endl;
		return 11 ;
	}
	if( dateMin<=0 ) {
		cout<<"Error. empty dateMin"<<endl;
		return 11 ;
	}
	if( dateMax<=0 ) {
		cout<<"Error. empty dateMax"<<endl;
		return 11 ;
	}
	if( blueMin > blueMax ) {
		cout<<"Error. blueMin greater than blueMax"<<endl;
		return 11 ;
	}
	if( outExtent_ullr.size()!=4 ) {
		cout<<"Error. invalid outExtent_ullr"<<endl;
		return 11 ;
	}
	if( outfilename =="") {
		cout<<"Error. empty outfilename."<<endl;
		return 11 ;
	}

	cout<<"---------"<<endl;
	cout<<"indir "<<indir <<endl;
	cout<<"prefix "<<fnprefix <<endl;
	cout<<"tail "<<fntail <<endl;
	cout<<"datepos "<<dpos <<endl;
	cout<<"datelen "<<dlen <<endl;
	cout<<"datemin "<<dateMin <<endl;
	cout<<"datemax "<<dateMax <<endl;
	cout<<"blueIndex0B "<<blueIndex0B <<endl;
	cout<<"vmin "<<blueMin <<endl;
	cout<<"vmax "<<blueMax <<endl;
	cout<<"fillvalue "<<fillvalue <<endl;
	cout<<"out ullr "<<outExtent_ullr[0]<<" "<<outExtent_ullr[1]<<" "<<outExtent_ullr[2]<<" "<<outExtent_ullr[3]<<endl;
	cout<<"outfilename "<<outfilename <<endl;
	cout<<"---------"<<endl;

	vector<string> infiles ;
	wDirTools wdir ;
	cout<<"Info: scan files..."<<endl;
	bool scandirok = wdir.getDirFiles( indir , fnprefix, fntail , infiles) ;
	if( scandirok == false )
	{
		//scan dir failed.
		cout<<"scan indir failed."<<endl;
		return 22 ;
	}
	cout<<"Info: scan files "<<infiles.size()<<endl;
	if( infiles.size()==0 ) {
		cout<<"Error. empty infiles"<<endl;
		return 22 ;
	}

	vector<string> usedFiles ;
	vector<int>    dateArray ;
	for(int iif = 0 ; iif < infiles.size() ;++iif) {
		string fullname = infiles[iif] ;
		string bname = wFilenameUtil::baseName(fullname);
		int datei = atof( bname.substr(dpos, dlen).c_str() ) ;
		if( datei >= dateMin && datei<=dateMax ) {
			usedFiles.push_back( fullname ) ;
			dateArray.push_back( datei) ;
		}
	}
	cout<<"Info: date filtered files "<<usedFiles.size()<<endl;
	if( usedFiles.size()==0 ) {
		cout<<"Error. empty filtered files."<<endl;
		return 22 ;
	}

	int inputBandCount = 0;
	double resolution = 0 ;
	string inputSrs ;
	GDALDataType dataType ;
	{//first input raster
		wGdalRaster* ptr = wGdalRasterFactory::OpenFile( usedFiles[0] ) ;
		if( ptr==0 ) {
			cout<<"Error. failed to open "<<usedFiles[0]<<endl;
			return 22 ;
		}
		inputBandCount = ptr->getNBand() ;
		const double* trans = ptr->getTrans() ;
		resolution = trans[1] ;
		inputSrs = string(ptr->getProj() );
		dataType = ptr->getDataType() ;
		delete ptr ;
	}
	double outTrans[6] = { outExtent_ullr[0], resolution, 0.0 , outExtent_ullr[1], 0.0 , -resolution } ;

	int outxsize = ceil((outExtent_ullr[2] - outExtent_ullr[0])/resolution) ;
	int outysize = ceil((outExtent_ullr[1] - outExtent_ullr[3])/resolution) ;

	cout<<"outsize "<<outxsize<<" "<<outysize<<endl;
	if( outxsize<=0 || outysize<=0 ) {
		cout<<"Error. out size invalid."<<endl;
		return 22 ;
	}
	int outbandCount = inputBandCount+1 ;
	const int tagBandIndex0B = inputBandCount;
	wGdalRaster* ptrOut = wGdalRasterFactory::Create(outxsize,outysize,outbandCount,dataType) ;
	if( ptrOut==0 ) {
		cout<<"Error. create output failed."<<endl;
		return 23 ;
	}
	for(int iy = 0 ; iy < outysize; ++ iy ) {
		for(int ix = 0 ; ix < outxsize; ++ ix ) {
			for(int iband=0;iband<inputBandCount;++iband) {
				ptrOut->setValued(ix,iy,iband, fillvalue) ;
			}
		}
	}

	for(int iin = 0 ; iin < usedFiles.size() ; ++ iin ) {
		cout<<"using " <<iin<<"/"<< usedFiles.size() <<" " <<usedFiles[iin]<<endl;
		wGdalRaster* inptr =  wGdalRasterFactory::OpenFile( usedFiles[iin] ) ;
		if( inptr == 0 ) {
			cout<<"Error. failed to open "<<usedFiles[iin]<<endl;
			delete ptrOut ;
			return 24 ;
		}
		int outTag = dateArray[iin]%10000 ;
		const double* trans = inptr->getTrans() ;
		int inxsize = inptr->getXSize() ;
		int inysize = inptr->getYSize() ;
		int nbands = inptr->getNBand() ;
		assert( nbands == inputBandCount) ;
		for(int iy = 0; iy < inysize; ++ iy ) {
			double incoory = trans[3] + trans[5] * iy + trans[5]/2.0 ;
			int outiy = (incoory - outTrans[3])/outTrans[5] ;
			if( outiy>=0 && outiy < outysize ){
				for(int ix=0;ix<inxsize;++ix) {
					double incoorx = trans[0] + trans[1] * ix + trans[1]/2.0 ;
					int outix = (incoorx - outTrans[0])/outTrans[1] ;
					if( outix >=0 && outix < outxsize ) {
						double oldbluevalue = ptrOut->getValued(outix,outiy,blueIndex0B) ;
						 
						double newbluevalue = inptr->getValued(ix,iy,blueIndex0B) ;
						if( newbluevalue>=blueMin && newbluevalue<=blueMax ) {
							if( newbluevalue < oldbluevalue || oldbluevalue==fillvalue ) {
								for(int iband=0;iband<inputBandCount;++iband) {
									ptrOut->setValued(outix,outiy,iband, inptr->getValued(ix,iy,iband)) ;
								}
								ptrOut->setValued(outix,outiy,tagBandIndex0B,outTag) ;
							}
						}
						
					}
				}
			}
		}
		delete inptr ; inptr = 0 ;
	}
	cout<<"saving "<<endl;
	ptrOut->copyTrans(outTrans);
	ptrOut->copyProj( inputSrs.c_str() ) ;
	bool saveok = ptrOut->save(outfilename) ; 
	if( saveok==false ) {
		delete ptrOut;
		ptrOut=0;
		cout<<"Error. failed to write "<<outfilename<<endl;
		return 30 ;
	}

	delete ptrOut;
	cout<<"done"<<endl;
	return 0 ;
}

